home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-19 / iritsm3s.zip / BSP_KNOT.C < prev    next >
C/C++ Source or Header  |  1991-10-09  |  17KB  |  496 lines
  1. /******************************************************************************
  2. * Bsp_knot.c - Various bspline routines to handle knot vectors.              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Aug. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #ifdef __MSDOS__
  8. #include <stdlib.h>
  9. #endif /* __MSDOS__ */
  10.  
  11. #include <ctype.h>
  12. #include <stdio.h>
  13. #include <string.h>
  14. #include "cagd_loc.h"
  15.  
  16. /******************************************************************************
  17. * Returns TRUE iff The given knot vector has open end conditions.          *
  18. ******************************************************************************/
  19. CagdBType BspKnotHasOpenEC(CagdRType *KnotVector, int Len, int Order)
  20. {
  21.     int i,
  22.         LastIndex = Len + Order - 1;
  23.  
  24.     for (i = 1; i < Order; i++)
  25.        if (!APX_EQ(KnotVector[0], KnotVector[i])) return FALSE;
  26.  
  27.     for (i = LastIndex - 1; i >= Len; i--)
  28.        if (!APX_EQ(KnotVector[LastIndex], KnotVector[i])) return FALSE;
  29.  
  30.     return TRUE;
  31. }
  32.  
  33. /******************************************************************************
  34. * Returns TRUE iff t is in the parametric domain as define by the knot vector *
  35. * KnotVector its length Len and the order Order.                  *
  36. ******************************************************************************/
  37. CagdBType BspKnotParamInDomain(CagdRType *KnotVector, int Len, int Order,
  38.                                    CagdRType t)
  39. {
  40.     CagdRType
  41.     r1 = KnotVector[Order - 1],
  42.     r2 = KnotVector[Len];
  43.  
  44.     return (r1 < t || APX_EQ(r1, t)) && (r2 > t || APX_EQ(r2, t));
  45. }
  46.  
  47. /******************************************************************************
  48. * Return the index of the last knot which is less than or equal to t in the   *
  49. * given knot vector KnotVector of length Len. APX_EQ is used for equality.    *
  50. * Parameter t is assumed to be in the parametric domain for the knot vector.  *
  51. ******************************************************************************/
  52. int BspKnotLastIndexLE(CagdRType *KnotVector, int Len, CagdRType t)
  53. {
  54.     int i;
  55.  
  56.     for (i = 0;
  57.      i < Len && (KnotVector[i] <= t || APX_EQ(KnotVector[i], t));
  58.      i++);
  59.  
  60.     return i - 1;
  61. }
  62.  
  63. /******************************************************************************
  64. * Return the index of the last knot which is less than t in the given knot    *
  65. * vector KnotVector of length Len.                          *
  66. * Parameter t is assumed to be in the parametric domain for the knot vector.  *
  67. ******************************************************************************/
  68. int BspKnotLastIndexL(CagdRType *KnotVector, int Len, CagdRType t)
  69. {
  70.     int i;
  71.  
  72.     for (i = 0;
  73.      i < Len && KnotVector[i] < t && !APX_EQ(KnotVector[i], t);
  74.      i++);
  75.  
  76.     return i - 1;
  77. }
  78.  
  79. /******************************************************************************
  80. * Return the index of the first knot which is greater than t in given knot    *
  81. * vector KnotVector of length Len.                          *
  82. * Parameter t is assumed to be in the parametric domain for the knot vector.  *
  83. ******************************************************************************/
  84. int BspKnotFirstIndexG(CagdRType *KnotVector, int Len, CagdRType t)
  85. {
  86.     int i;
  87.  
  88.     for (i = Len - 1;
  89.      i >= 0 && KnotVector[i] > t && !APX_EQ(KnotVector[i], t);
  90.      i--);
  91.  
  92.     return i + 1;
  93. }
  94.  
  95. /******************************************************************************
  96. * Return a uniform float knot vector for Len Control points and order Order.  *
  97. * If KnotVector is NULL it is being allocated dynamically.              *
  98. ******************************************************************************/
  99. CagdRType *BspKnotUniformFloat(int Len, int Order, CagdRType *KnotVector)
  100. {
  101.     int i;
  102.     CagdRType *KV;
  103.  
  104.     if (KnotVector == NULL)
  105.     KnotVector = KV = (CagdRType *) CagdMalloc(sizeof(CagdRType) *
  106.                                 (Len + Order));
  107.     else
  108.     KV = KnotVector;
  109.  
  110.     for (i = 0; i < Len + Order; i++) *KV++ = (CagdRType) i;
  111.  
  112.     return KnotVector;
  113. }
  114.  
  115. /******************************************************************************
  116. * Return a uniform open knot vector for Len Control points and order Order.   *
  117. * If KnotVector is NULL it is being allocated dynamically.              *
  118. ******************************************************************************/
  119. CagdRType *BspKnotUniformOpen(int Len, int Order, CagdRType *KnotVector)
  120. {
  121.     int i, j;
  122.     CagdRType *KV;
  123.  
  124.     if (KnotVector == NULL)
  125.     KnotVector = KV = (CagdRType *) CagdMalloc(sizeof(CagdRType) *
  126.                                 (Len + Order));
  127.     else
  128.     KV = KnotVector;
  129.  
  130.     for (i = 0; i < Order; i++) *KV++ = 0.0;
  131.     for (i = 1, j = Len - Order; i <= j; ) *KV++ = (CagdRType) i++;
  132.     for (j = 0; j < Order; j++) *KV++ = (CagdRType) i;
  133.  
  134.     return KnotVector;
  135. }
  136.  
  137. /******************************************************************************
  138. * Apply an affine transformation to the given knot vector. Note affine        *
  139. * transformation for the knot vector does not change the spline.          *
  140. * Knot vector is translated by Translate amount and scaled by Scale.
  141. ******************************************************************************/
  142. void BspKnotAffineTrans(CagdRType *KnotVector, int Len,
  143.                     CagdRType Translate, CagdRType Scale)
  144. {
  145.     int i;
  146.  
  147.     for (i = 0; i < Len; i++)
  148.     KnotVector[i] = (KnotVector[i] + Translate) * Scale;
  149. }
  150.  
  151. /******************************************************************************
  152. * Creates an identical copy of a given knot vector KnotVector of length Len.  *
  153. ******************************************************************************/
  154. CagdRType *BspKnotCopy(CagdRType *KnotVector, int Len)
  155. {
  156.     CagdRType
  157.     *NewKnotVector = (CagdRType *) CagdMalloc(sizeof(CagdRType) * Len);
  158.  
  159.     GEN_COPY(NewKnotVector, KnotVector, sizeof(CagdRType) * Len);
  160.  
  161.     return NewKnotVector;
  162. }
  163.  
  164. /******************************************************************************
  165. * Reverse a knot vector of length Len.                          *
  166. ******************************************************************************/
  167. CagdRType *BspKnotReverse(CagdRType *KnotVector, int Len)
  168. {
  169.     int i;
  170.     CagdRType
  171.     *NewKnotVector = (CagdRType *) CagdMalloc(sizeof(CagdRType) * Len),
  172.     t = KnotVector[Len - 1] + KnotVector[0];
  173.  
  174.     for (i = 0; i < Len; i++)
  175.     NewKnotVector[i] = t - KnotVector[Len - i - 1];
  176.  
  177.     return NewKnotVector;
  178. }
  179.  
  180. /******************************************************************************
  181. * Returns all the knots in KnotVector1 not in KnotVector2.              *
  182. * NewLen is set to new KnotVector length.                      *
  183. ******************************************************************************/
  184. CagdRType *BspKnotSubtrTwo(CagdRType *KnotVector1, int Len1,
  185.                CagdRType *KnotVector2, int Len2,
  186.                int *NewLen)
  187. {
  188.     int i = 0,
  189.     j = 0;
  190.     CagdRType
  191.     *NewKnotVector = (CagdRType *) CagdMalloc(sizeof(CagdRType) * Len1),
  192.     *t = NewKnotVector;
  193.  
  194.     *NewLen = 0;
  195.     while (i < Len1 && j < Len2) {
  196.     if (APX_EQ(KnotVector1[i], KnotVector2[j])) {
  197.         i++;
  198.         j++;
  199.     }
  200.     else if (KnotVector1[i] > KnotVector2[j]) {
  201.         j++;
  202.     }
  203.     else {
  204.         /* Current KnotVector1 value is less than current KnotVector2: */
  205.         *t++ = KnotVector1[i++];
  206.         (*NewLen)++;
  207.     }
  208.     }
  209.  
  210.     return NewKnotVector;
  211. }
  212.  
  213. /******************************************************************************
  214. * Merge two knot vector KnotVector1 and KnotVector2 of length Len1 and Len2   *
  215. * respectively into one. If Mult is not zero then knot multiplicity is tested *
  216. * not to be bigger than Mult value. NewLen is set to new KnotVector length.   *
  217. ******************************************************************************/
  218. CagdRType *BspKnotMergeTwo(CagdRType *KnotVector1, int Len1,
  219.                CagdRType *KnotVector2, int Len2,
  220.                int Mult, int *NewLen)
  221. {
  222.     int i = 0,
  223.     j = 0,
  224.     k = 0,
  225.         m = 0,
  226.     Len = Len1 + Len2;
  227.     CagdRType t,
  228.     *NewKnotVector = (CagdRType *) CagdMalloc(sizeof(CagdRType) * Len);
  229.  
  230.     if (Mult == 0) Mult = Len + 1;
  231.  
  232.     while (i < Len1 && j < Len2) {
  233.     if (KnotVector1[i] < KnotVector2[j]) {
  234.         /* Use value from KnotVector1: */
  235.         t = KnotVector1[i++];
  236.     }
  237.     else {
  238.         /* Use value from KnotVector2: */
  239.         t = KnotVector2[j++];
  240.     }
  241.  
  242.     if (k == 0 || k > 0 && !APX_EQ(NewKnotVector[k - 1], t)) {
  243.         NewKnotVector[k++] = t;
  244.         m = 1;
  245.     }
  246.     else if (m < Mult) {
  247.         NewKnotVector[k++] = t;
  248.         m++;
  249.     }
  250.     }
  251.  
  252.     /* It should be noted that k <= Len so there is a chance some of the new */
  253.     /* KnotVector space will not be used (it is not memory leak!). If you    */
  254.     /* really care that much - copy it to a new allocated vector of size k   */
  255.     /* and return it while freeing the original of size Len.             */
  256.     *NewLen = k;
  257.  
  258.     return NewKnotVector;
  259. }
  260.  
  261. /******************************************************************************
  262. * Creates a new vector with the given KnotVector Node values.              *
  263. * The nodes are the approximated parametric value associated with the each    *
  264. * control point. Therefore for a knot vector of length Len and order Order    *
  265. * there are Len - Order control points and therefore nodes.              *
  266. * Nodes are defined as (k = Order - 1 or degree):                  *
  267. *                                              *
  268. *      i+k                                      *
  269. *       -                 First Node N(i = 0)              *
  270. *       \                                      *
  271. *          / KnotVector(j)        Last Node  N(i = Len - k - 2)          *
  272. *       -                                      *
  273. *        j=i+1                                      *
  274. * N(i) = -----------------                              *
  275. *            k                                  *
  276. ******************************************************************************/
  277. CagdRType *BspKnotNodes(CagdRType *KnotVector, int Len, int Order)
  278. {
  279.     int i, j,
  280.     k = Order - 1,
  281.     NodeLen = Len - Order;
  282.     CagdRType
  283.     *NodeVector = (CagdRType *) CagdMalloc(sizeof(CagdRType) * NodeLen);
  284.  
  285.     for (i = 0; i < NodeLen; i++) {
  286.     NodeVector[i] = 0.0;
  287.     for (j = i + 1; j <= i + k; j++) NodeVector[i] += KnotVector[j];
  288.     NodeVector[i] /= k;
  289.     }
  290.  
  291.     return NodeVector;
  292. }
  293.  
  294. /******************************************************************************
  295. * Creates a new vector with t inserted as a new knot. Attempt is made to make *
  296. * sure t is in the knot vector domain.                          *
  297. * No test is made for the current multiplicity for knot t in KnotVector.      *
  298. ******************************************************************************/
  299. CagdRType *BspKnotInsertOne(CagdRType *KnotVector, int Order, int Len,
  300.                                 CagdRType t)
  301. {
  302.     int Mult = BspKnotFindMult(KnotVector, Order, Len, t) + 1;
  303.  
  304.     return BspKnotInsertMult(KnotVector, Order, &Len, t, Mult);
  305. }
  306.  
  307. /******************************************************************************
  308. * Inserts Mult knots with value t into the knot vector KnotVector.          *
  309. * Attempt is made to make sure t in knot vector domain.                  *
  310. * If a knot equal to t (up to APX_EQ) already exists with multiplicity i only *
  311. * Mult - i knot are been inserted into the new knot vector.              *
  312. * Len is updated to the resulting knot vector.                      *
  313. * Note it is possible to DELETE a knot using this routine by specifying       *
  314. * multiplicity less then current multiplicity!                      *
  315. ******************************************************************************/
  316. CagdRType *BspKnotInsertMult(CagdRType *KnotVector, int Order, int *Len,
  317.                                 CagdRType t, int Mult)
  318. {
  319.     int i, CurrentMult, NewLen, FirstIndex;
  320.     CagdRType *NewKnotVector;
  321.  
  322.     if (!BspKnotParamInDomain(KnotVector, *Len, Order, t))
  323.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  324.  
  325.     CurrentMult = BspKnotFindMult(KnotVector, Order, *Len, t);
  326.     NewLen = *Len + Mult - CurrentMult;
  327.     NewKnotVector = (CagdRType *) CagdMalloc(sizeof(CagdRType) *
  328.                                 (NewLen + Order));
  329.     FirstIndex = BspKnotLastIndexL(KnotVector, *Len + Order, t) + 1;
  330.  
  331.     /* Copy all the knot before the knot t. */
  332.     GEN_COPY(NewKnotVector, KnotVector, sizeof(CagdRType) * FirstIndex);
  333.  
  334.     /* Insert Mult knot of value t. */
  335.     for (i = FirstIndex; i < FirstIndex + Mult; i++)
  336.     NewKnotVector[i] = t;
  337.  
  338.     /* And copy the second part. */
  339.     GEN_COPY(&NewKnotVector[FirstIndex + Mult],
  340.          &KnotVector[FirstIndex + CurrentMult],
  341.          sizeof(CagdRType) * (*Len + Order - FirstIndex - Mult + 1));
  342.  
  343.     *Len = NewLen;
  344.     return NewKnotVector;
  345. }
  346.  
  347. /******************************************************************************
  348. * Returns the multiplicity of knot t in knot vector KnotVector, zero if none. *
  349. ******************************************************************************/
  350. int BspKnotFindMult(CagdRType *KnotVector, int Order, int Len, CagdRType t)
  351. {
  352.     int LastIndex, FirstIndex;
  353.  
  354.     if (!BspKnotParamInDomain(KnotVector, Len, Order, t))
  355.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  356.  
  357.     FirstIndex = BspKnotLastIndexL(KnotVector, Len, t) + 1;
  358.  
  359.     for (LastIndex = FirstIndex;
  360.      LastIndex < Len && APX_EQ(KnotVector[LastIndex], t);
  361.      LastIndex++);
  362.  
  363.     return LastIndex - FirstIndex;
  364. }
  365.  
  366. /******************************************************************************
  367. * Scans the given knot vector to a potential C1 discontinuity.              *
  368. *   Returns TRUE if found one and set t to its parameter value.              *
  369. *   Assumes knot vector has open end condition.                      *
  370. ******************************************************************************/
  371. CagdBType BspKnotC1Discont(CagdRType *KnotVector, int Order, int Length,
  372.                                 CagdRType *t)
  373. {
  374.     int i, Count;
  375.     CagdRType
  376.     LastT = KnotVector[0] - 1.0;
  377.  
  378.     for (i = Order, Count = 0; i < Length; i++) {
  379.     if (APX_EQ(LastT, KnotVector[i]))
  380.         Count++;
  381.     else {
  382.         Count = 1;
  383.         LastT = KnotVector[i];
  384.     }
  385.  
  386.     if (Count >= Order - 1) {
  387.         *t = LastT;
  388.         return TRUE;
  389.     }
  390.     }
  391.  
  392.     return FALSE;
  393. }
  394.  
  395. /******************************************************************************
  396. * Scans the given knot vector to aall potential C1 discontinuity. Returns     *
  397. * a vector holding the parameter values of C1 discontinuities, NULL of non    *
  398. * found.                                      *
  399. *   Sets n to length of returned vector.                      *
  400. *   Assumes knot vector has open end condition.                      *
  401. ******************************************************************************/
  402. CagdRType *BspKnotAllC1Discont(CagdRType *KnotVector, int Order, int Length,
  403.                                     int *n)
  404. {
  405.     int i, Count,
  406.     C1DiscontCount = 0;
  407.     CagdRType
  408.     LastT = KnotVector[0] - 1.0,
  409.     *C1Disconts = (CagdRType *) CagdMalloc(sizeof(CagdRType) * Length);
  410.  
  411.     for (i = Order, Count = 0; i < Length; i++) {
  412.     if (APX_EQ(LastT, KnotVector[i]))
  413.         Count++;
  414.     else {
  415.         Count = 1;
  416.         LastT = KnotVector[i];
  417.     }
  418.  
  419.     if (Count >= Order - 1) {
  420.         C1Disconts[C1DiscontCount++] = LastT;
  421.         Count = 0;
  422.     }
  423.     }
  424.  
  425.     if ((*n = C1DiscontCount) > 0) {
  426.     return C1Disconts;
  427.     }
  428.     else {
  429.     CagdFree((VoidPtr) C1Disconts);
  430.     return NULL;
  431.     }
  432. }
  433.  
  434. /*****************************************************************************
  435. * Routine to determine where to sample along the provided paramtric domain,  *
  436. * given the C1 discontinuities along it.                     *
  437. *   Returns a vector of length NumSamples.                     *
  438. *   If C1Disconts != NULL (NumC1Disconts > 0) this vector is being freed.    *
  439. *****************************************************************************/
  440. CagdRType *BspKnotParamValues(CagdRType PMin, CagdRType PMax, int NumSamples,
  441.                   CagdRType *C1Disconts, int NumC1Disconts)
  442. {
  443.     int i, CrntIndex, NextIndex, CrntDiscont;
  444.     CagdRType Step,
  445.     *Samples = (CagdRType *) CagdMalloc(sizeof(CagdRType) * NumSamples);
  446.  
  447.     Samples[0] = PMin;
  448.     if (NumSamples <= 1) return Samples;
  449.     Samples[NumSamples - 1] = PMax;
  450.     if (NumSamples <= 2) return Samples;
  451.  
  452.     if (NumC1Disconts > NumSamples - 2) {
  453.         /* More C1 discontinuities than we are sampling. Grab a subset of    */
  454.     /* The discontinuities as the sampling set.                 */
  455.         Step = ((CagdRType) (NumC1Disconts - 1)) / (NumSamples - 2) - EPSILON;
  456.     for (i = 0; i < NumSamples - 2; i++)
  457.             Samples[i + 1] = C1Disconts[(int) (i * Step)];
  458.     }
  459.     else {
  460.     /* More samples than C1 discontinuites. Uniformly distribute the C1  */
  461.     /* discontinuites between the samples and linearly interpolate the   */
  462.     /* sample values in between.                         */
  463.         Step = ((CagdRType) (NumC1Disconts + 1)) / (NumSamples - 2);
  464.     CrntIndex = CrntDiscont = 0;
  465.  
  466.     while (CrntDiscont < NumC1Disconts) {
  467.         NextIndex = (int) ((CrntDiscont + 1) / Step);
  468.         Samples[NextIndex] = C1Disconts[CrntDiscont++];
  469.  
  470.         for (i = CrntIndex + 1; i < NextIndex; i++) {
  471.         CagdRType
  472.             t = (i - CrntIndex) / ((CagdRType) (NextIndex - CrntIndex)),
  473.             t1 = 1.0 - t;
  474.  
  475.         Samples[i] = Samples[CrntIndex] * t1 + Samples[NextIndex] * t;
  476.         }
  477.  
  478.         CrntIndex = NextIndex;
  479.     }
  480.  
  481.     /* Interpolate the last interval: */
  482.         for (i = CrntIndex + 1; i < NumSamples - 1; i++) {
  483.             CagdRType
  484.             t = (i - CrntIndex) / ((CagdRType) (NumSamples - 1 - CrntIndex)),
  485.             t1 = 1.0 - t;
  486.  
  487.             Samples[i] = Samples[CrntIndex] * t1 + Samples[NumSamples - 1] * t;
  488.         }
  489.     }
  490.  
  491.     if (C1Disconts != NULL)
  492.     CagdFree((VoidPtr) C1Disconts);
  493.  
  494.     return Samples;
  495. }
  496.